CESM2-Large Ensemble Reproduction of Kay et al. 2015¶
Introduction¶
This Jupyter Notebook demonstrates how one might use the NCAR Community Earth System Model v2 (CESM2) Large Ensemble (LENS) data hosted on AWS S3. The notebook shows how to reproduce figures 2 and 4 from the Kay et al. (2015) paper describing the CESM LENS dataset (doi:10.1175/BAMS-D-13-00255.1), with the LENS2 dataset.
This resource is intended to be helpful for people not familiar with elements of the Pangeo framework including Jupyter Notebooks, Xarray, and Zarr data format, or with the original paper, so it includes additional explanation.
Imports¶
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
# Silence dask.distributed logs
import logging
logger = logging.getLogger("distributed.utils_perf")
logger.setLevel(logging.ERROR)
import intake
import numpy as np
import pandas as pd
import xarray as xr
import hvplot.pandas, hvplot.xarray
import holoviews as hv
from ncar_jobqueue import NCARCluster
from distributed import Client
hv.extension('bokeh')
Spin up a Cluster¶
cluster = NCARCluster(memory="50GB", walltime='2:00:00', cores=4, processes=1, resource_spec='select=1:ncpus=1:mem=50GB')
cluster.scale(20)
client = Client(cluster)
client
Client
Client-7310d807-2633-11ec-848e-3cecef1b11fa
| Connection method: Cluster object | Cluster type: dask_jobqueue.PBSCluster |
| Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/mgrover/proxy/37785/status |
Cluster Info
PBSCluster
f562885f
| Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/mgrover/proxy/37785/status | Workers: 0 |
| Total threads: 0 | Total memory: 0 B |
Scheduler Info
Scheduler
Scheduler-2d673e6f-43b9-4cce-95ac-2d66638dc1b5
| Comm: tcp://10.12.206.54:33624 | Workers: 0 |
| Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/mgrover/proxy/37785/status | Total threads: 0 |
| Started: Just now | Total memory: 0 B |
Workers
Read in Data¶
Use the Intake-ESM Catalog to Access the Data¶
catalog = intake.open_esm_datastore(
'../data/glade-cesm2-le.json'
)
catalog
aws-cesm2-le catalog with 23 dataset(s) from 205 asset(s):
| unique | |
|---|---|
| component | 3 |
| spatial_domain | 2 |
| frequency | 3 |
| path | 196 |
| experiment | 2 |
| forcing_variant | 2 |
| variable | 35 |
| long_name | 35 |
| start_time | 4 |
| end_time | 4 |
Subset for Daily Temperature Data (TREFHT)¶
We use Intake-ESM here to query for our dataset, focusing the smoothed biomass burning (smbb) experiment!
catalog_subset = catalog.search(variable='TREFHT', frequency='daily')
catalog_subset
aws-cesm2-le catalog with 4 dataset(s) from 4 asset(s):
| unique | |
|---|---|
| component | 1 |
| spatial_domain | 1 |
| frequency | 1 |
| path | 4 |
| experiment | 2 |
| forcing_variant | 2 |
| variable | 1 |
| long_name | 1 |
| start_time | 2 |
| end_time | 2 |
catalog_subset.df
| component | spatial_domain | frequency | path | experiment | forcing_variant | variable | long_name | start_time | end_time | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | historical | cmip6 | TREFHT | Reference height temperature | 1850-01-01 12:00:00 | 2015-01-01 12:00:00 |
| 1 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | historical | smbb | TREFHT | Reference height temperature | 1850-01-01 12:00:00 | 2015-01-01 12:00:00 |
| 2 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | ssp370 | cmip6 | TREFHT | Reference height temperature | 2015-01-01 12:00:00 | 2101-01-01 12:00:00 |
| 3 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | ssp370 | smbb | TREFHT | Reference height temperature | 2015-01-01 12:00:00 | 2101-01-01 12:00:00 |
Taking a look at the dataframe, we see there are two files - one for the historical run, and the other a future scenario (ssp370)
catalog_subset.df
| component | spatial_domain | frequency | path | experiment | forcing_variant | variable | long_name | start_time | end_time | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | historical | cmip6 | TREFHT | Reference height temperature | 1850-01-01 12:00:00 | 2015-01-01 12:00:00 |
| 1 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | historical | smbb | TREFHT | Reference height temperature | 1850-01-01 12:00:00 | 2015-01-01 12:00:00 |
| 2 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | ssp370 | cmip6 | TREFHT | Reference height temperature | 2015-01-01 12:00:00 | 2101-01-01 12:00:00 |
| 3 | atm | global | daily | /glade/scratch/mgrover/data/lens2-aws/atm/dail... | ssp370 | smbb | TREFHT | Reference height temperature | 2015-01-01 12:00:00 | 2101-01-01 12:00:00 |
Load in our datasets using .to_dataset_dict()¶
We use .to_dataset_dict() to return a dictionary of datasets, which we can now use for the analysis
dsets = catalog_subset.to_dataset_dict()
--> The keys in the returned dictionary of datasets are constructed as follows:
'component.experiment.frequency.forcing_variant'
Here are the keys for each dataset - you will notice they follow the groupby attributes, with component.experiment.frequency.forcing_variant, one for the future and the other historical here
dsets.keys()
dict_keys(['atm.ssp370.daily.smbb', 'atm.historical.daily.smbb', 'atm.historical.daily.cmip6', 'atm.ssp370.daily.cmip6'])
Let’s load these into separate datasets, then merge these together!
historical_smbb = dsets['atm.historical.daily.smbb']
future_smbb = dsets['atm.ssp370.daily.smbb']
historical_cmip6 = dsets['atm.historical.daily.cmip6']
future_cmip6 = dsets['atm.ssp370.daily.cmip6']
merge_ds_smbb = xr.concat([historical_smbb, future_smbb], dim='time')
merge_ds_smbb = merge_ds_smbb.dropna(dim='member_id')
merge_ds_cmip6= xr.concat([historical_cmip6, future_cmip6], dim='time')
merge_ds_cmip6 = merge_ds_cmip6.dropna(dim='member_id')
merge_ds_smbb['time'] = merge_ds_smbb.indexes['time'].to_datetimeindex()
merge_ds_cmip6['time'] = merge_ds_cmip6.indexes['time'].to_datetimeindex()
Finally, we can load in just the variable we are interested in, Reference height temperature (TREFHT). We will separate this into two arrays:
t_smbb- the entire period of data for the smoothed biomass burning experiment (1850 to 2100)t_cmip6- the entire period of data for the cmip6 biomass burning experiment (1850 to 2100)t_ref- the reference period used in the Kay et al. 2015 paper (1961 to 1990)
t_smbb = merge_ds_smbb.TREFHT
t_cmip6 = merge_ds_cmip6.TREFHT
t_ref = t_cmip6.sel(time=slice('1961', '1990'))
Read in the Grid Data¶
We also have Zarr stores of grid data, such as the area of each grid cell. We will need to follow a similar process here, querying for our experiment and extracting the dataset
grid_subset = catalog.search(component='atm', frequency='static', experiment='historical', forcing_variant='smbb')
_, grid = grid_subset.to_dataset_dict(aggregate=False, zarr_kwargs={"consolidated": True}).popitem()
grid
--> The keys in the returned dictionary of datasets are constructed as follows:
'component.spatial_domain.frequency.path.experiment.forcing_variant.variable.long_name.start_time.end_time'
<xarray.Dataset>
Dimensions: (lat: 192, lon: 288, ilev: 31, lev: 30, bnds: 2, slat: 191, slon: 288)
Coordinates: (12/20)
P0 float64 ...
area (lat, lon) float32 dask.array<chunksize=(192, 288), meta=np.ndarray>
gw (lat) float64 dask.array<chunksize=(192,), meta=np.ndarray>
hyai (ilev) float64 dask.array<chunksize=(31,), meta=np.ndarray>
hyam (lev) float64 dask.array<chunksize=(30,), meta=np.ndarray>
hybi (ilev) float64 dask.array<chunksize=(31,), meta=np.ndarray>
... ...
ntrm int32 ...
ntrn int32 ...
* slat (slat) float64 -89.53 -88.59 -87.64 -86.7 ... 87.64 88.59 89.53
* slon (slon) float64 -0.625 0.625 1.875 3.125 ... 355.6 356.9 358.1
w_stag (slat) float64 dask.array<chunksize=(191,), meta=np.ndarray>
wnummax (lat) int32 dask.array<chunksize=(192,), meta=np.ndarray>
Dimensions without coordinates: bnds
Data variables:
*empty*
Attributes:
intake_esm_varname: None
intake_esm_dataset_key: atm.global.static./glade/scratch/mgrover/data/le...- lat: 192
- lon: 288
- ilev: 31
- lev: 30
- bnds: 2
- slat: 191
- slon: 288
- P0()float64...
- long_name :
- reference pressure
- units :
- Pa
array(100000.)
- area(lat, lon)float32dask.array<chunksize=(192, 288), meta=np.ndarray>
- long_name :
- Grid-Cell Area
- standard_name :
- cell_area
- units :
- m2
Array Chunk Bytes 216.00 kiB 216.00 kiB Shape (192, 288) (192, 288) Count 2 Tasks 1 Chunks Type float32 numpy.ndarray - gw(lat)float64dask.array<chunksize=(192,), meta=np.ndarray>
- long_name :
- gauss weights
Array Chunk Bytes 1.50 kiB 1.50 kiB Shape (192,) (192,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - hyai(ilev)float64dask.array<chunksize=(31,), meta=np.ndarray>
- long_name :
- hybrid A coefficient at layer interfaces
Array Chunk Bytes 248 B 248 B Shape (31,) (31,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - hyam(lev)float64dask.array<chunksize=(30,), meta=np.ndarray>
- long_name :
- hybrid A coefficient at layer midpoints
Array Chunk Bytes 240 B 240 B Shape (30,) (30,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - hybi(ilev)float64dask.array<chunksize=(31,), meta=np.ndarray>
- long_name :
- hybrid B coefficient at layer interfaces
Array Chunk Bytes 248 B 248 B Shape (31,) (31,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - hybm(lev)float64dask.array<chunksize=(30,), meta=np.ndarray>
- long_name :
- hybrid B coefficient at layer midpoints
Array Chunk Bytes 240 B 240 B Shape (30,) (30,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - ilev(ilev)float642.255 5.032 10.16 ... 985.1 1e+03
- formula_terms :
- a: hyai b: hybi p0: P0 ps: PS
- long_name :
- hybrid level at interfaces (1000*(A+B))
- positive :
- down
- standard_name :
- atmosphere_hybrid_sigma_pressure_coordinate
- units :
- level
array([ 2.25524 , 5.031692, 10.157947, 18.555317, 30.669123, 45.867477, 63.323483, 80.701418, 94.941042, 111.693211, 131.401271, 154.586807, 181.863353, 213.952821, 251.704417, 296.117216, 348.366588, 409.835219, 482.149929, 567.224421, 652.332969, 730.445892, 796.363071, 845.353667, 873.715866, 900.324631, 924.964462, 947.432335, 967.538625, 985.11219 , 1000. ]) - lat(lat)float64-90.0 -89.06 -88.12 ... 89.06 90.0
- axis :
- Y
- bounds :
- lat_bnds
- long_name :
- latitude
- standard_name :
- latitude
- units :
- degrees_north
array([-90. , -89.057594, -88.115181, -87.172775, -86.23037 , -85.287956, -84.345551, -83.403145, -82.460732, -81.518326, -80.575912, -79.633507, -78.691101, -77.748688, -76.806282, -75.863876, -74.921463, -73.979057, -73.036652, -72.094238, -71.151833, -70.209427, -69.267014, -68.324608, -67.382202, -66.439789, -65.497383, -64.554977, -63.612564, -62.670158, -61.727749, -60.785339, -59.842934, -58.900524, -57.958115, -57.015705, -56.073299, -55.13089 , -54.18848 , -53.246075, -52.303665, -51.361256, -50.41885 , -49.47644 , -48.534031, -47.591621, -46.649216, -45.706806, -44.764397, -43.821991, -42.879581, -41.937172, -40.994766, -40.052357, -39.109947, -38.167538, -37.225132, -36.282722, -35.340313, -34.397907, -33.455498, -32.513088, -31.570681, -30.628273, -29.685863, -28.743456, -27.801046, -26.858639, -25.916231, -24.973822, -24.031414, -23.089005, -22.146597, -21.204189, -20.26178 , -19.319372, -18.376963, -17.434555, -16.492147, -15.549738, -14.607329, -13.664922, -12.722513, -11.780105, -10.837696, -9.895288, -8.95288 , -8.010471, -7.068063, -6.125654, -5.183246, -4.240838, -3.298429, -2.356021, -1.413613, -0.471204, 0.471204, 1.413613, 2.356021, 3.298429, 4.240838, 5.183246, 6.125654, 7.068063, 8.010471, 8.95288 , 9.895288, 10.837696, 11.780105, 12.722513, 13.664922, 14.607329, 15.549738, 16.492147, 17.434555, 18.376963, 19.319372, 20.26178 , 21.204189, 22.146597, 23.089005, 24.031414, 24.973822, 25.916231, 26.858639, 27.801046, 28.743456, 29.685863, 30.628273, 31.570681, 32.513088, 33.455498, 34.397907, 35.340313, 36.282722, 37.225132, 38.167538, 39.109947, 40.052357, 40.994766, 41.937172, 42.879581, 43.821991, 44.764397, 45.706806, 46.649216, 47.591621, 48.534031, 49.47644 , 50.41885 , 51.361256, 52.303665, 53.246075, 54.18848 , 55.13089 , 56.073299, 57.015705, 57.958115, 58.900524, 59.842934, 60.785339, 61.727749, 62.670158, 63.612564, 64.554977, 65.497383, 66.439789, 67.382202, 68.324608, 69.267014, 70.209427, 71.151833, 72.094238, 73.036652, 73.979057, 74.921463, 75.863876, 76.806282, 77.748688, 78.691101, 79.633507, 80.575912, 81.518326, 82.460732, 83.403145, 84.345551, 85.287956, 86.23037 , 87.172775, 88.115181, 89.057594, 90. ]) - lat_bnds(lat, bnds)float64dask.array<chunksize=(192, 2), meta=np.ndarray>
Array Chunk Bytes 3.00 kiB 3.00 kiB Shape (192, 2) (192, 2) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - lev(lev)float643.643 7.595 14.36 ... 976.3 992.6
- formula_terms :
- a: hyam b: hybm p0: P0 ps: PS
- long_name :
- hybrid level at midpoints (1000*(A+B))
- positive :
- down
- standard_name :
- atmosphere_hybrid_sigma_pressure_coordinate
- units :
- level
array([ 3.643466, 7.59482 , 14.356632, 24.61222 , 38.2683 , 54.59548 , 72.012451, 87.82123 , 103.317127, 121.547241, 142.994039, 168.22508 , 197.908087, 232.828619, 273.910817, 322.241902, 379.100904, 445.992574, 524.687175, 609.778695, 691.38943 , 763.404481, 820.858369, 859.534767, 887.020249, 912.644547, 936.198398, 957.48548 , 976.325407, 992.556095]) - lon(lon)float640.0 1.25 2.5 ... 356.2 357.5 358.8
- axis :
- X
- bounds :
- lon_bnds
- long_name :
- longitude
- standard_name :
- longitude
- units :
- degrees_east
array([ 0. , 1.25, 2.5 , ..., 356.25, 357.5 , 358.75])
- lon_bnds(lon, bnds)float64dask.array<chunksize=(288, 2), meta=np.ndarray>
Array Chunk Bytes 4.50 kiB 4.50 kiB Shape (288, 2) (288, 2) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - ntrk()int32...
- long_name :
- spectral truncation parameter K
array(1, dtype=int32)
- ntrm()int32...
- long_name :
- spectral truncation parameter M
array(1, dtype=int32)
- ntrn()int32...
- long_name :
- spectral truncation parameter N
array(1, dtype=int32)
- slat(slat)float64-89.53 -88.59 ... 88.59 89.53
- long_name :
- staggered latitude
- units :
- degrees_north
array([-89.528796, -88.586387, -87.643979, -86.701571, -85.759162, -84.816754, -83.874346, -82.931937, -81.989529, -81.04712 , -80.104712, -79.162304, -78.219895, -77.277487, -76.335079, -75.39267 , -74.450262, -73.507853, -72.565445, -71.623037, -70.680628, -69.73822 , -68.795812, -67.853403, -66.910995, -65.968586, -65.026178, -64.08377 , -63.141361, -62.198953, -61.256545, -60.314136, -59.371728, -58.429319, -57.486911, -56.544503, -55.602094, -54.659686, -53.717277, -52.774869, -51.832461, -50.890052, -49.947644, -49.005236, -48.062827, -47.120419, -46.17801 , -45.235602, -44.293194, -43.350785, -42.408377, -41.465969, -40.52356 , -39.581152, -38.638743, -37.696335, -36.753927, -35.811518, -34.86911 , -33.926702, -32.984293, -32.041885, -31.099476, -30.157068, -29.21466 , -28.272251, -27.329843, -26.387435, -25.445026, -24.502618, -23.560209, -22.617801, -21.675393, -20.732984, -19.790576, -18.848168, -17.905759, -16.963351, -16.020942, -15.078534, -14.136126, -13.193717, -12.251309, -11.308901, -10.366492, -9.424084, -8.481675, -7.539267, -6.596859, -5.65445 , -4.712042, -3.769634, -2.827225, -1.884817, -0.942408, 0. , 0.942408, 1.884817, 2.827225, 3.769634, 4.712042, 5.65445 , 6.596859, 7.539267, 8.481675, 9.424084, 10.366492, 11.308901, 12.251309, 13.193717, 14.136126, 15.078534, 16.020942, 16.963351, 17.905759, 18.848168, 19.790576, 20.732984, 21.675393, 22.617801, 23.560209, 24.502618, 25.445026, 26.387435, 27.329843, 28.272251, 29.21466 , 30.157068, 31.099476, 32.041885, 32.984293, 33.926702, 34.86911 , 35.811518, 36.753927, 37.696335, 38.638743, 39.581152, 40.52356 , 41.465969, 42.408377, 43.350785, 44.293194, 45.235602, 46.17801 , 47.120419, 48.062827, 49.005236, 49.947644, 50.890052, 51.832461, 52.774869, 53.717277, 54.659686, 55.602094, 56.544503, 57.486911, 58.429319, 59.371728, 60.314136, 61.256545, 62.198953, 63.141361, 64.08377 , 65.026178, 65.968586, 66.910995, 67.853403, 68.795812, 69.73822 , 70.680628, 71.623037, 72.565445, 73.507853, 74.450262, 75.39267 , 76.335079, 77.277487, 78.219895, 79.162304, 80.104712, 81.04712 , 81.989529, 82.931937, 83.874346, 84.816754, 85.759162, 86.701571, 87.643979, 88.586387, 89.528796]) - slon(slon)float64-0.625 0.625 1.875 ... 356.9 358.1
- long_name :
- staggered longitude
- units :
- degrees_east
array([ -0.625, 0.625, 1.875, ..., 355.625, 356.875, 358.125])
- w_stag(slat)float64dask.array<chunksize=(191,), meta=np.ndarray>
- long_name :
- staggered latitude weights
Array Chunk Bytes 1.49 kiB 1.49 kiB Shape (191,) (191,) Count 2 Tasks 1 Chunks Type float64 numpy.ndarray - wnummax(lat)int32dask.array<chunksize=(192,), meta=np.ndarray>
- long_name :
- cutoff Fourier wavenumber
Array Chunk Bytes 768 B 768 B Shape (192,) (192,) Count 2 Tasks 1 Chunks Type int32 numpy.ndarray
- intake_esm_varname :
- None
- intake_esm_dataset_key :
- atm.global.static./glade/scratch/mgrover/data/lens2-aws/atm/static/grid.zarr.historical.smbb.nan.nan.nan.nan
Extract the Grid Area and Compute the Total Area¶
Here, we extract the grid area and compute the total area across all grid cells - this will help when computing the weights…
ds_grid = xr.open_dataset('/glade/campaign/cgd/cesm/CESM2-LE/timeseries/atm/proc/tseries/month_1/AREA/b.e21.BHISTcmip6.f09_g17.LE2-1001.001.cam.h0.AREA.185001-185912.nc')
cell_area = ds_grid.isel(time=0).AREA.load()
total_area = cell_area.sum()
cell_area
<xarray.DataArray 'AREA' (lat: 192, lon: 288)>
array([[2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07,
2.9948368e+07, 2.9948368e+07],
[2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08,
2.3957478e+08, 2.3957478e+08],
[4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08,
4.7908477e+08, 4.7908477e+08],
...,
[4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08,
4.7908477e+08, 4.7908477e+08],
[2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08,
2.3957478e+08, 2.3957478e+08],
[2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07,
2.9948368e+07, 2.9948368e+07]], dtype=float32)
Coordinates:
* lat (lat) float64 -90.0 -89.06 -88.12 -87.17 ... 87.17 88.12 89.06 90.0
* lon (lon) float64 0.0 1.25 2.5 3.75 5.0 ... 355.0 356.2 357.5 358.8
time object 1850-02-01 00:00:00
Attributes:
units: m2
long_name: area of grid box
cell_methods: time: mean- lat: 192
- lon: 288
- 2.995e+07 2.995e+07 2.995e+07 ... 2.995e+07 2.995e+07 2.995e+07
array([[2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07, 2.9948368e+07, 2.9948368e+07], [2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08, 2.3957478e+08, 2.3957478e+08], [4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08, 4.7908477e+08, 4.7908477e+08], ..., [4.7908477e+08, 4.7908477e+08, 4.7908477e+08, ..., 4.7908477e+08, 4.7908477e+08, 4.7908477e+08], [2.3957478e+08, 2.3957478e+08, 2.3957478e+08, ..., 2.3957478e+08, 2.3957478e+08, 2.3957478e+08], [2.9948368e+07, 2.9948368e+07, 2.9948368e+07, ..., 2.9948368e+07, 2.9948368e+07, 2.9948368e+07]], dtype=float32) - lat(lat)float64-90.0 -89.06 -88.12 ... 89.06 90.0
- long_name :
- latitude
- units :
- degrees_north
array([-90. , -89.057592, -88.115183, -87.172775, -86.230366, -85.287958, -84.34555 , -83.403141, -82.460733, -81.518325, -80.575916, -79.633508, -78.691099, -77.748691, -76.806283, -75.863874, -74.921466, -73.979058, -73.036649, -72.094241, -71.151832, -70.209424, -69.267016, -68.324607, -67.382199, -66.439791, -65.497382, -64.554974, -63.612565, -62.670157, -61.727749, -60.78534 , -59.842932, -58.900524, -57.958115, -57.015707, -56.073298, -55.13089 , -54.188482, -53.246073, -52.303665, -51.361257, -50.418848, -49.47644 , -48.534031, -47.591623, -46.649215, -45.706806, -44.764398, -43.82199 , -42.879581, -41.937173, -40.994764, -40.052356, -39.109948, -38.167539, -37.225131, -36.282723, -35.340314, -34.397906, -33.455497, -32.513089, -31.570681, -30.628272, -29.685864, -28.743455, -27.801047, -26.858639, -25.91623 , -24.973822, -24.031414, -23.089005, -22.146597, -21.204188, -20.26178 , -19.319372, -18.376963, -17.434555, -16.492147, -15.549738, -14.60733 , -13.664921, -12.722513, -11.780105, -10.837696, -9.895288, -8.95288 , -8.010471, -7.068063, -6.125654, -5.183246, -4.240838, -3.298429, -2.356021, -1.413613, -0.471204, 0.471204, 1.413613, 2.356021, 3.298429, 4.240838, 5.183246, 6.125654, 7.068063, 8.010471, 8.95288 , 9.895288, 10.837696, 11.780105, 12.722513, 13.664921, 14.60733 , 15.549738, 16.492147, 17.434555, 18.376963, 19.319372, 20.26178 , 21.204188, 22.146597, 23.089005, 24.031414, 24.973822, 25.91623 , 26.858639, 27.801047, 28.743455, 29.685864, 30.628272, 31.570681, 32.513089, 33.455497, 34.397906, 35.340314, 36.282723, 37.225131, 38.167539, 39.109948, 40.052356, 40.994764, 41.937173, 42.879581, 43.82199 , 44.764398, 45.706806, 46.649215, 47.591623, 48.534031, 49.47644 , 50.418848, 51.361257, 52.303665, 53.246073, 54.188482, 55.13089 , 56.073298, 57.015707, 57.958115, 58.900524, 59.842932, 60.78534 , 61.727749, 62.670157, 63.612565, 64.554974, 65.497382, 66.439791, 67.382199, 68.324607, 69.267016, 70.209424, 71.151832, 72.094241, 73.036649, 73.979058, 74.921466, 75.863874, 76.806283, 77.748691, 78.691099, 79.633508, 80.575916, 81.518325, 82.460733, 83.403141, 84.34555 , 85.287958, 86.230366, 87.172775, 88.115183, 89.057592, 90. ]) - lon(lon)float640.0 1.25 2.5 ... 356.2 357.5 358.8
- long_name :
- longitude
- units :
- degrees_east
array([ 0. , 1.25, 2.5 , ..., 356.25, 357.5 , 358.75])
- time()object1850-02-01 00:00:00
- long_name :
- time
- bounds :
- time_bnds
array(cftime.DatetimeNoLeap(1850, 2, 1, 0, 0, 0, 0, has_year_zero=True), dtype=object)
- units :
- m2
- long_name :
- area of grid box
- cell_methods :
- time: mean
Compute the Weighted Annual Averages¶
Here, we utilize the resample(time="AS") method which does an annual resampling based on start of calendar year.
Note we also weight by the cell area, multiplying each value by the corresponding cell area, summing over the grid (lat, lon), then dividing by the total area
Setup the Computation - Below, we lazily prepare the calculation!¶
You’ll notice how quickly this cell runs - which is suspicious 👀 - we didn’t actually do any computation, but rather prepared the calculation
%%time
t_ref_ts = (
(t_ref.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon'))
).mean(dim=("time", "member_id"))
t_smbb_ts = (
(t_smbb.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon')))
t_cmip6_ts = (
(t_cmip6.resample(time="AS").mean("time")).weighted(cell_area).mean(('lat', 'lon')))
CPU times: user 1.95 s, sys: 14.7 ms, total: 1.96 s
Wall time: 1.95 s
t_ref_ts.values
array(287.47504, dtype=float32)
Compute the Averages¶
At this point, we actually compute the averages - using %%time to help us measure how long it takes to compute…
%%time
t_ref_mean = t_ref_ts.compute()
CPU times: user 15.3 s, sys: 514 ms, total: 15.8 s
Wall time: 20.6 s
%%time
t_smbb_mean = t_smbb_ts.compute()
CPU times: user 1min 42s, sys: 3.94 s, total: 1min 46s
Wall time: 2min 3s
%%time
t_cmip6_mean = t_cmip6_ts.compute()
CPU times: user 2min 11s, sys: 4.99 s, total: 2min 16s
Wall time: 2min 42s
Convert to Dataframes¶
Our output data format is still an xarray.DataArray which isn’t neccessarily ideal… we could convert this to a dataframe which would be easier to work with!
We can convert our xarray.DataArray to a pandas dataframe using the following:
t_smbb_ts_df = t_smbb_mean.to_series().unstack().T
t_cmip6_ts_df = t_cmip6_mean.to_series().unstack().T
anomaly_smbb = (t_smbb_ts_df - t_ref_mean.data)
anomaly_cmip6 = (t_cmip6_ts_df - t_ref_mean.data)
Grab some Observational Data¶
In this case, we use the HADCRUT4 dataset
# Observational time series data for comparison with ensemble average
obsDataURL = "https://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/cru/hadcrut4/air.mon.anom.median.nc"
ds = xr.open_dataset(obsDataURL).load()
ds
<xarray.Dataset>
Dimensions: (lat: 36, lon: 72, time: 2060, nbnds: 2)
Coordinates:
* lat (lat) float32 87.5 82.5 77.5 72.5 ... -72.5 -77.5 -82.5 -87.5
* lon (lon) float32 -177.5 -172.5 -167.5 -162.5 ... 167.5 172.5 177.5
* time (time) datetime64[ns] 1850-01-01 1850-02-01 ... 2021-08-01
Dimensions without coordinates: nbnds
Data variables:
time_bnds (time, nbnds) datetime64[ns] 1850-01-01 1850-01-31 ... 2021-08-31
air (time, lat, lon) float32 nan nan nan nan nan ... nan nan nan nan
Attributes:
platform: Surface
title: HADCRUT4 Combined Air Temperature/SST An...
history: Originally created at NOAA/ESRL PSD by C...
Conventions: CF-1.0
Comment: This dataset supersedes V3
Source: Obtained from http://hadobs.metoffice.co...
version: 4.2.0
dataset_title: HadCRUT4
References: https://www.psl.noaa.gov/data/gridded/da...
DODS_EXTRA.Unlimited_Dimension: time- lat: 36
- lon: 72
- time: 2060
- nbnds: 2
- lat(lat)float3287.5 82.5 77.5 ... -82.5 -87.5
- long_name :
- Latitude
- units :
- degrees_north
- actual_range :
- [ 87.5 -87.5]
- axis :
- Y
- coordinate_defines :
- point
- standard_name :
- latitude
- _ChunkSizes :
- 36
array([ 87.5, 82.5, 77.5, 72.5, 67.5, 62.5, 57.5, 52.5, 47.5, 42.5, 37.5, 32.5, 27.5, 22.5, 17.5, 12.5, 7.5, 2.5, -2.5, -7.5, -12.5, -17.5, -22.5, -27.5, -32.5, -37.5, -42.5, -47.5, -52.5, -57.5, -62.5, -67.5, -72.5, -77.5, -82.5, -87.5], dtype=float32) - lon(lon)float32-177.5 -172.5 ... 172.5 177.5
- long_name :
- Longitude
- units :
- degrees_east
- actual_range :
- [-177.5 177.5]
- axis :
- X
- coordinate_defines :
- point
- standard_name :
- longitude
- _ChunkSizes :
- 72
array([-177.5, -172.5, -167.5, -162.5, -157.5, -152.5, -147.5, -142.5, -137.5, -132.5, -127.5, -122.5, -117.5, -112.5, -107.5, -102.5, -97.5, -92.5, -87.5, -82.5, -77.5, -72.5, -67.5, -62.5, -57.5, -52.5, -47.5, -42.5, -37.5, -32.5, -27.5, -22.5, -17.5, -12.5, -7.5, -2.5, 2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 37.5, 42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5, 87.5, 92.5, 97.5, 102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5, 152.5, 157.5, 162.5, 167.5, 172.5, 177.5], dtype=float32) - time(time)datetime64[ns]1850-01-01 ... 2021-08-01
- long_name :
- Time
- delta_t :
- 0000-01-00 00:00:00
- avg_period :
- 0000-01-00 00:00:00
- standard_name :
- time
- axis :
- T
- coordinate_defines :
- start
- bounds :
- time_bnds
- actual_range :
- [18262. 80931.]
- _ChunkSizes :
- 1
array(['1850-01-01T00:00:00.000000000', '1850-02-01T00:00:00.000000000', '1850-03-01T00:00:00.000000000', ..., '2021-06-01T00:00:00.000000000', '2021-07-01T00:00:00.000000000', '2021-08-01T00:00:00.000000000'], dtype='datetime64[ns]')
- time_bnds(time, nbnds)datetime64[ns]1850-01-01 ... 2021-08-31
- long_name :
- Time Boundaries
- _ChunkSizes :
- [1 2]
array([['1850-01-01T00:00:00.000000000', '1850-01-31T00:00:00.000000000'], ['1850-02-01T00:00:00.000000000', '1850-02-28T00:00:00.000000000'], ['1850-03-01T00:00:00.000000000', '1850-03-31T00:00:00.000000000'], ..., ['2021-06-01T00:00:00.000000000', '2021-06-30T00:00:00.000000000'], ['2021-07-01T00:00:00.000000000', '2021-07-31T00:00:00.000000000'], ['2021-08-01T00:00:00.000000000', '2021-08-31T00:00:00.000000000']], dtype='datetime64[ns]') - air(time, lat, lon)float32nan nan nan nan ... nan nan nan nan
- dataset :
- HADCRUT4
- var_desc :
- Air Temperature
- level_desc :
- Surface
- statistic :
- Anomaly
- parent_stat :
- Observation
- valid_range :
- [-40. 40.]
- units :
- degC
- long_name :
- HADCRUT4: Median Surface Air Temperature Monthly Median Anomaly from 100 ensemble members
- precision :
- 2
- cell_methods :
- time: anomaly (monthly from values)
- standard_name :
- air_temperature_anomaly
- actual_range :
- [-9.9692100e+36 2.0620308e+01]
- _ChunkSizes :
- [ 1 36 72]
array([[[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ..., [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], ... [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]], [[ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, 0.3049999 , ..., nan, nan, nan], [ nan, nan, 1.8495195 , ..., nan, nan, nan], ..., [ nan, nan, nan, ..., -2.7343245 , nan, nan], [ nan, nan, nan, ..., nan, nan, nan], [ nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
- platform :
- Surface
- title :
- HADCRUT4 Combined Air Temperature/SST Anomaly
- history :
- Originally created at NOAA/ESRL PSD by CAS 04/2012 from files obtained at the Hadley Center
- Conventions :
- CF-1.0
- Comment :
- This dataset supersedes V3
- Source :
- Obtained from http://hadobs.metoffice.com/ Data is a collaboration of CRU and the Hadley Center
- version :
- 4.2.0
- dataset_title :
- HadCRUT4
- References :
- https://www.psl.noaa.gov/data/gridded/data.hadcrut4.html
- DODS_EXTRA.Unlimited_Dimension :
- time
Compute the Weighted Temporal Mean from Seasons¶
def weighted_temporal_mean(ds):
"""
weight by days in each month
"""
time_bound_diff = ds.time_bnds.diff(dim="nbnds")[:, 0]
wgts = time_bound_diff.groupby("time.year") / time_bound_diff.groupby(
"time.year"
).sum(xr.ALL_DIMS)
np.testing.assert_allclose(wgts.groupby("time.year").sum(xr.ALL_DIMS), 1.0)
obs = ds["air"]
cond = obs.isnull()
ones = xr.where(cond, 0.0, 1.0)
obs_sum = (obs * wgts).resample(time="AS").sum(dim="time")
ones_out = (ones * wgts).resample(time="AS").sum(dim="time")
obs_s = (obs_sum / ones_out).mean(("lat", "lon")).to_series()
return obs_s
obs_s = weighted_temporal_mean(ds)
Convert the dataset into a dataframe¶
obs_df = pd.DataFrame(obs_s).rename(columns={0:'value'})
obs_df
| value | |
|---|---|
| time | |
| 1850-01-01 | -0.338822 |
| 1851-01-01 | -0.245482 |
| 1852-01-01 | -0.291014 |
| 1853-01-01 | -0.342457 |
| 1854-01-01 | -0.276820 |
| ... | ... |
| 2017-01-01 | 0.777498 |
| 2018-01-01 | 0.641953 |
| 2019-01-01 | 0.809306 |
| 2020-01-01 | 0.865248 |
| 2021-01-01 | 0.683412 |
172 rows × 1 columns
Plot the Output¶
In this this case, we use hvPlot to plot the output!
smbb_plot = (anomaly_smbb.hvplot.scatter('time', height=300, width=500, color='lightgrey', legend=False, xlabel='Year', ylabel='Global Mean \n Temperature Anomaly (K)', hover=False) *\
anomaly_smbb.mean(1).hvplot.line('time', color='k', line_width=3, label='ensemble mean') *\
obs_df.hvplot.line('time', color='red', label='observations')).opts(title='Smoothed Biomass Burning')
cmip6_plot = anomaly_cmip6.hvplot.scatter('time', height=300, width=500, color='lightgrey', legend=False, xlabel='Year', ylabel='Global Mean \n Temperature Anomaly (K)', hover=False) *\
anomaly_cmip6.mean(1).hvplot.line('time', color='k', line_width=3, label='ensemble mean') *\
obs_df.hvplot.line('time', color='red', label='observations').opts(title='CMIP6 Biomass Burning')
smbb_plot + cmip6_plot
Spin Down the Cluster¶
After we are done, we can spin down our cluster
cluster.close()
client.close()